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Abstract. The Diffusion Monte Carlo method is devoted to the computation of electronic 
ground-state energies of molecules. In this paper, we focus on implementations of this method 
which consist in exploring the configuration space with a fixed number of random walkers 
evolving according to a Stochastic Differential Equation discretized in time. We allow stochastic 
reconfigurations of the walkers to reduce the discrepancy between the weights that they carry. 
On a simple one-dimensional example, we prove the convergence of the method for a fixed 
number of reconfigurations when the number of walkers tends to +00 while the timestep tends 
to 0. We confirm our theoretical rates of convergence by numerical experiments. Various 
resampling algorithms are investigated, both theoretically and numerically. 
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Introduction 

The computation of electronic structures of atoms, molecules and solids is a central problem in chemistry 
and physics. We focus here on electronic ground state calculations where the objective is the computation 
of the lowest eigenvalue (the so-called ground-state energy) Eq of a self-adjoint Hamiltonian H = — ^A+V 
with domain D-h(H) on a Hilbcrt space H C L 2 (M. 3N ) where N is the number of electrons (see [4] for a 
general introduction) : 

E = inf{(V>,ffV>, i> 6 D n (H), U\\ = 1}, (1) 

where (•, •) denotes the duality bracket on L 2 (M. 3N ) and || • || the L 2 (R 3Ar )-norm. For simplicity, we 
omit the spin variables. The function V describes the interaction between the electrons, and between 
the electrons and the nuclei, which are supposed to be fixed point-like particles. The functions ip are 
square integrable, their normalized square modulus \ip\ 2 being interpreted as the probability density of the 
particles positions in space, and they satisfy an antisymmetry condition with respect to the numbering 
of the electrons, due to the fermionic nature of the electrons (Pauli principle): H = Ai=i £ 2 (R 3 )- We 
suppose that the potential V is such that Eq is an isolated eigenvalue of H (see [3] for sufficient conditions), 
and we denote by tpo a normalized eigenfunction associated with Eq. 
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Due to the high dimensionality of the problem, stochastic methods are particularly well suited to com- 
pute Eq. The first particle approximation scheme of such spectral quantities was introduced in [12] for 
finite state space models. Convergence analysis for such interacting particle systems (both continuous or 
discrete in time) first appeared in [7-10]. The Diffusion Monte Carlo (DMC) method is widely used in 
chemistry (see [2,17]), but has been only recently considered from a mathematical viewpoint (see [3,14]). 
This method gives an estimate of Eq in terms of the long-time limit of the expectation of a functional 
of a drift-diffusion process with a source term. It requires an importance sampling function ipj which 
approximates the ground-state i/jq of H. Let us define the drift function b = Vln\ipi\, the so-called local 

energy El = and the DMC energy: 

ipi 

E (E L {X t ) cxp (- /„* E L {X s )ds^ } 7 
£ dmcW = 7 7 rx -, (2) 

where the 3TV-dimensional process X t satisfies the stochastic differential equation: 

(3) 



| X t = X Q + J b(X s )ds + dW t 
\ X Q ^\^\ 2 (x)dx. 



The stochastic process (W t )t>o is a standard 3TV-dimcnsional Brownian motion. One can then show 
that (see [3]) 

lim EuMc(t) = £dmc,o, (4) 

t— >oo 

where 

Bdmc,o = m{{(i>,Hi>), t/j G D n (H), = 1, ip = on ^(O)}. (5) 

We have proved in [3] that EuMCfi > Eq, with equality if and only if the nodal surfaces of ipi coincide 
with those of a ground state i/>o of H. In other words, if there exists a ground state "0o such that 
ij)J {$) = ijjQ (0), then lim^oo £7dmc(£) = E a . The error \E a — -Edmc,o| is related to the so-called 
fixed-node approximation, which is well known by practitioners of the field (sec [4]). 

In this paper, we complement the theoretical results obtained in [3] with a numerical analysis in a simple 
case. In practice, the longtime limit -Edmco i n Q is approximated by taking the value of -Edmc at a 
(large) time T > 0. Then -Edmc(^) is approximated by using a discretization in time of the stochastic 
differential equation ([3]) and of the integral in the exponential factor in ^ , and an approximation of the 
expectation values in ^ by an empirical mean over a large number N of trajectories. These trajectories 
(X l )i<a<N, also called walkers in the physical literature or particles in the mathematical literature, satisfy 
a discretized version of ([3]), and interact at times nAt for n € {1, . . . ,v — 1} where At = T/v for i/£M' 
through a stochastic reconfiguration step aimed at reducing the discrepancy between their exponential 
weights. We thus obtain an interacting particle system. The number of reconfiguration steps is v — \. The 
stochastic differential equation ((3]) is discretized with a possibly smaller timestep 8t = At/n = T/(vn) 
with k£N*. The total number of steps for the discretization of is then K = vk. 

In the following, we consider the following adapted version of the DMC scheme with a fixed number of 
walkers (see [2]): 

• Initialization of an ensemble of TV walkers (-X"oAt) i.i.d. according to \ip i\ 2 {x) dx . 

V / l<j<N 

• Iterations in time: let us be given the particle positions (_X^ At ) at time nAt, for 

V nii V l<j<N 

n G {0, . . . , v — 1}. The new particle positions at time (n + l)At are obtained in two steps: 

(1) Walkers displacement: for all 1 < j < TV, the successive positions (x J nAt+st , . . . , X 3 nAt+KS 

over the time interval (nAt, (n+ l)At) are obtained by an appropriate discretization of 
In the field of interacting particles system for Feynman-Kac formulae (sec [7.9]). this step is 
called the mutation step. 
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(2) Stochastic reconfiguration: The new positions^ (xi n+1 ^ At \ which will be used as 



>+ 1 ) A V i<j<;v 

the initial particle positions on the time interval ((n + l)Ai, (n + 2)Ai) are obtained from 
independent sampling of the measure 



E^Li ex P f-^ELi E d x l&t+k8t) 



In words, the new particle positions ( Xj , nAt ) arc randomly chosen among the 



final particle positions ( X nAt+KSt J , each of them being weighted with the coefficient 



cxp y— <5<Efc=i E L{X J nAt+kSt )J (accordingly to the exponential factor in |J3j)). In the field 
of interacting particles system for Feynman-Kac formulae, this step is called the selection 
step. 

An estimate of i?DMc(^n+i) is then given by: 

1 N 

EBMG(tn+l) ~ E L ( X l n+1)At ) • ( 7 ) 

J = l 

There are other possible estimations of £x)Mc(in+i) • ln [2], the authors propose to use Cesaro or weighted 
Cesaro means of the expression (|7|). In Section [TJ we will use the following expression: 

EjLi ^£(^«Ai+«5t) ex p(- ft Efc=l £; i( X nAt+fc5t)) 
^DMCl^ri+l) — ~ 7 " : \ j (8) 

At+kSt) I 

in an intermediate step to prove the convergence result. 

We would like to mention that a continuous in time version of the DMC scheme with stochastic reconfig- 
uration has been proposed in [14]. The author analyzes the longtime behavior of the interacting particle 
system and proves in particular a uniform in time control of the variance of the estimated energy. 

The DMC algorithm presented above is prototypical. Many refinements are used in practice. For exam- 
ple, an acception-rejection step is generally used in the walkers displacement step (see [13]). This will not 
be discussed here. Likewise, the selection step can be done in many ways (see [5,6] for general algorithms, 
and [2, 15, 17] for algorithms used in the context of DMC computations). In this paper, we restrict our- 
selves to resampling methods with a fixed number of particles, and such that the weights of the particles 
after resampling are equal to 1. Then, the basic consistency requirement of the selection step is that, 
conditionally on the former positions I X J nAt , kSt J , the i-th particle X l nAt+KSt is replicated 

V / l<j<N,l<k<K 

Np l n times in mean, where p l n = exp (-^ELi E h{ X \At+tet)) / E^Li ex P (-^ELi E d xJ n At+ksS) 
denotes the (normalized) weight of the i-th particle. There arc of course many ways to satisfy this re- 
quirement. We presented above the so-called multinomial resampling method. We will also discuss below 
residual resampling (also called stochastic remainder resampling) , stratified resampling and systematic re- 
sampling, which may also be used for DMC computations. Let us briefly describe these three resampling 
methods. Residual resampling consists in reproducing |_iVp^J times the i-th particle, and then completing 
the set of particles by using multinomial resampling to draw the N R = N — Ejli [Np l n \ remaining parti- 
cles, the i-th particle being assigned the weight p R > 1 = {Np l n }/N R . Here and in the following, [x\ and {x} 
respectively denote the integer and the fractional part of x £ K. In the stratified resampling method, the 



^With a slight abuse of notation and though nAt + re<5i = (n + l)At, we distinguish between the particle positions 
-^nAt+nSt a * * ne ent ^ °f the walkers displacement on time interval (nAt, (n+l)At), and the new particle positions ^^ n+1 ^ t 
obtained after the reconfiguration step, and which are used as the initial position for the next walkers displacement on time 
interval ((ra + l)At, (ra + 2) At). We will use a more precise notation for the analysis of the numerical scheme in Section [T] 
but this is not required at this stage. 
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interval (0, 1) is divided into N intervals ((i—\)/N,i/N) (1 < i < N), N random variables are then drawn 
independently and uniformly in each interval, and the new particle positions are then obtained by the 
inversion method: X\ n+l)At = J2f=i P l n <(i-u} l )/ N <'Ei=i p l n } X li5t+nSv whcrc U n arc Li - d - random 

variables uniformly distributed over [0, 1]. Here and in the following, we use the convention Ym=i ' = 0- 
Systematic resampling consists in replicating the z'-th particle NJ2l=i Pn + U n — Pn + U n 

time^, where (U n ) n >i are independent random variables uniformly distributed in [0, 1]. Notice that sys 
tematic resampling can be seen as the stratified resampling method, with Jj\ = . . . = f/„ = U n . Contrary 
to the three other resampling methods, after a systematic resampling step, the new particle positions are 
not independent, conditionally on the former positions. This makes systematic resampling much more 
difficult to study mathematically. To our knowledge, its convergence even in a discrete time setting is 
still an open question. We will therefore restrict ourselves to a numerical study of its performance. 

Notice that practitioners often use branching algorithms with an evolving number of walkers during the 
computation (see [13,17]): the particles with low local energy are replicated and the particles with high 
local energy are killed, without keeping the total number of particles constant. This may lead to a smaller 
Monte Carlo error (fourth contribution to the error in the classification just below). 



1 N 

We can distinguish between four sources of errors in the approximation of Eq by — El [ 

.7=1 



'3 



(1) the error due to the fixed node approximation \E — -EdmcoIi 

(2) the error due to finite time approximation of the limit: limt_>oo £"dmc(^) — £dmc(T)i 

(3) the error due to the time discretization of the stochastic differential equation ([3]) and of the 
integral in the exponential factor in EuMc(t) (see ([2])), 

(4) the error introduced by the interacting particle system, due to the approximation of the expec- 
tation value in ([2]) by an empirical mean. 

The error (1) due to the fixed node approximation has been analyzed theoretically in [3]. 

Concerning the error (2) due to finite time approximation of the limit, the rate of convergence in time 
is typically exponential. Indeed if H admits a spectral gap (namely if the distance between Eq and the 
remaining of the spectrum of H is strictly positive), and if ipi is such that (?/>/, Hip r ) < inf a css (H), then 
one can show that the operator H with domain Du{H) n {ip, ip = on ■0J 1 (O)} (whose lowest eigenvalue 
is £dmc,o, see ((5])) also admits a spectral gap 7 > 0. Then, by standard spectral decomposition methods, 
we have: 

< |£dmcW - £dmc,o| < Ceaq>(-7t). 

Our aim in this paper is to provide some theoretical and numerical results related to the errors (3) 
and (4), in the framework of a simple one-dimensional case. We therefore consider in the following that 
the final time of simulation T is fixed and we analyze the error introduced by the numerical scheme on 
the estimate of Eumc(T). Our convergence result is of the form: 



E 



1 N 

Eumc(T) - E L (X 



i=i 



' VK&t 



<C {T)5t+^j^, (9) 



where G(T) (resp. C(T,v)) denotes a constant which only depends on T (rcsp. on T and v) (see 
Theorem |4] and Corollary [T3l below) . 

Let us now present the toy model we consider in the following. We consider the Hamiltonian 

H = ~-^+V, with V = ^-x 2 + 6x\ (10) 
2 ax z 2 



2 The consistency of this resampling method follows from the following easy computation 

E([x + U\) = [x\P(U < l-{x}) + {[x\ +1)F(U > 1 - {x}) = [x\{l-{x}) + ([x\ +l){x}=x. 
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where u, 9 > are two constants. The ground state energy Eq is defined by {T]), with 

H = {tp G L 2 (R), ip(x) = -ip(-x)} . (11) 

We restrict the functional spaces to odd functions in order to mimic the antisymmetry constraint on tp for 
fermionic systems. The importance sampling ipi is chosen to be the ground state of Hq = — \-^z + \x 2 
on H: 

1 /4 

ip^x) = \F2w (-J xe~^ x \ (12) 

It is associated with the energy |w: H^ipi = ^Loipi- The drift function b and the local energy El are 
then defined by: 

b(x) = tl(x) = ~- uox, and E L (x) = V(x) - I ^-(x) = + 6x A . (13) 
ipi x 2-0/2 

Thus, using equation @, the DMC energy is: 

E DMC (t) = -u + 9^- j^— » (14) 

E(e X p(-9f*Xfds)) 



where 



X t =X + J^ (J^-uX s ^Jds + Wt, (15) 

with (Wt)t>o a Brownian motion independent from the initial variable Xq which is distributed according 
to the invariant measure 2ipj(x)l{ x>0 ydx. We recall that due to the explosive part in the drift function b, 
the stochastic process cannot cross 0, which is the zero point of 0/ (see [3]): ¥(3t > 0, X t = 0) = 0. This 
explains why the restriction of ip] to is indeed an invariant measure for (|15[) . For 6 > 0, the longtime 
limit i?DMC,o °f Ei)Mc{t) is not analytically known, but can be very accurately computed by a spectral 
method (see Section 12.11) . Let us finally make precise that for the numerical analysis, we use a special 
feature of our simple model, namely the fact that for s < t, it is possible to simulate the conditional law 
of X t given X s (see Appendix). The time discretization error is thus only related to the discretization 
of the integral in the exponential factor in the DMC energy ([2]). We however indicate some possible 
ways to prove {]]) with a convenient time discretization of the SDE (see Equation (|17p . Remark [3] and 
Proposition [T4|) . 

Though our model is one-dimensional (and therefore still far from the real problem |l])), it contains one of 
the main difficulties related to the approximation of the ground state energy for fermionic systems, namely 
the explosive behavior of the drift in the stochastic differential equation. However, two characteristics of 
practical problems are missing in the toy model considered here. First, since we consider a one-particle 
model, we do not treat difficulties related to singularities of the drift and of the local energy at points 
where two particles (either two electrons or one electron and one nucleus) coincide. Second, the local 
energy El generally explodes at the nodes of the trial wave function, and this is not the case on the 
simple example we study since the trial wave function is closely related to the exact ground state. For 
an adaptation of the DMC algorithm to take care of these singularities, we refer to [17]. Despite the 
simplicity of the model studied in this paper, we think that the convergence results we obtain and the 
mathematical tools we use are prototypical for generalization to more complicated systems. 

Compared to previous mathematical analysis of convergence for interacting particle systems with stochas- 
tic reconfiguration [7-10, 14], our study concentrates on the limit St — > and N — > oo for a fixed time T, 
and on the influence of the time discretization error in the estimate (j9|), where the test function El is 
unbounded. It is actually important in our analysis that this unbounded function El also appears in the 
weights of the particles, since it allows for specific estimates (see Lemmas 1^1 and [Til below) . 

The paper is organized as follows. In Section[TJ we prove the convergence result, by adapting the methods 
of [7, 9] to analyze the dependence of the error on St. We then check the optimality of this theoretical 
result by numerical experiments in Section [21 where we also analyze numerically the dependence of the 
results on various numerical parameters, including the number {y — 1) of reconfiguration steps. From 
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these numerical experiments, we propose a simple heuristic method to choose the optimal number of 
reconfiguration steps. 

Notation: For any set of random variables (1^)^/, we denote by cr((li)i e j) the sigma-ficld generated by 
these random variables. The parameters u> and 9 are fixed positive constants. By convention, any sum 
from one to zero is equal to zero: Efc=i ' = 0- Likewise, the subset {1, 2, . . . , 0} of N is by convention the 
empty set. For any real x, \_x\ and {x} respectively denote the integer and the fractional part of x. 



1. Numerical Analysis in a Simple Case 

We perform the numerical analysis in two steps: time discretization and then particle approximation. 



1.1. Time discretization 

We recall that T > denotes the final simulation time, and that St = ^ is the smallest time-step. Since 
Y t = X? is a square root process solving dY t = (3 — 2u)Yt)dt + 2^/YfdWt, it is possible to simulate the 
increments Y^+\)st — Ykst and therefore X(k+i)5t — XkSt (see Appendix or [11] p. 120). We can thus 
simulate exactly in law the vector (Xq, Xgt, ■ ■ ■ , Xxst)- That is why we are first going to study the error 
related to the time discretization of the integral which appears in the exponential factors in (I14|) . 



Let us define the corresponding approximation of Edmc{T)'- 
Eg MC (T) = ^-, 7 v - xV^ = £ w + *— ^ ~" ' " . (16) 



E[E L (X T )exp[-5tY:ti E L(X kSt ))) 3 E (Xj, exp (—6 St J2k=i X kSt) 



E (ocp \-8t Ef =1 E L (X k5t )j J 2 E (exp [ -9St ££ =1 X^ 

Proposition 1. 

VJT € N*, \E UMC (T) - Eg uc (T)\ < C T 6t 



Proof : Using Holder inequality, we have: 

(x£exp (-0 Jo Xfds 



\E BMC (T) -Eg MC (T)\ <— _ y/E&*) 



(exp(-^Ef=i^)) V E(exp(-0/ o T ^ 
exp (-9 J X*dsj - exp (-96t £ X 4 kSt \ ^ 



The conclusion is now a consequence of Lemma [2] and the fact that the function x e K + — > e 6x 
Lipschitz continuous with constant 9. 



Lemma 2. For any K G N* ; 



e ( (j* xjds -stJ2 **«J I < c&ey* + t) 



2\ 



where St = . 
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Proof of Lemma H : By Ito's formula, dXf = (10A 2 - 4uXf)dt + iXfdW t . With the integration by 
parts formula, one deduces that for any k € {1, . . . , A'}, 

kSt pkSt 



(X* st - Xf)ds = (s-(k- l)St) ((10X 2 S - 4 W A s 4 )ds + 4A S W S ) . 

(k-l)St J(k-l)5t 

Therefore denoting r s = L^J*^ the discretization time just before s, one obtains 

StY^Xtst - / Xfds = (s- r s )(10X s 2 - AuX A s )ds + (s - r s )4X*dW s . 
k=1 Jo Jo Jo 

Hence 

E U**S X fc«-/ T ^ ds ) j <2^ T ( S -r s ) 2 E(T(10X s 2 -4^ s 4 ) 2 + 16X s 6 ))d S . 

Since Ao is distributed according to the invariant measure 2ip](x)l{ x> o}dx, so is X s . As a consequence, 
for any pgN, E(AP) does not depend on s and is finite and the conclusion follows readily. 



In realistic situations, exact simulation of the increments Xi k+1 -jg t — X k g t is not possible and one has 
to resort to discretization schemes. The singularity of the drift coefficient prevents the process A t from 
crossing the nodal surfaces of the importance sampling function The standard explicit Eulcr scheme 
does not preserve this property at the discretized level. For that purpose, we suggest to use the following 
explicit scheme proposed by [1] 



Ao — A , 

Vfc e N, X {k+m 



X k5t {l-u8t) 



A%i 
1-uSt 



1/2 



25t 



with AW, 



k+l 



(fc+l)<St 



W, 



kSt ■ 



(17) 

Because of the singularity at the origin of the drift coefficient in (|15| . we have not been able so far to 
prove the following weak error bound (see Remark [3] below) : 



E /(A|)exp [-6 / A> 



K 



/(X|)exp l-OSt^X 



kSi 



k=l 



Such a bound is expected according to [16] and would imply that 



< C T St for f(x) = 1 and x 4 . 

(18) 



E 



Edmc(T) — 



(e l {X t ) exp (-St J2k =1 E L (X kSt ) 



cxp ( StJ2 k=1 E L (X kSt ) 



< c T st. 



(19) 



Remark 3. We would like to sketch a possible way to prove (| 18|) . Because the square root in (|17[) makes 
expansions with respect to St and AWk+i complicated, it is easier to work with Y t — A 2 and Y k $t = X kSt 
which satisfy 



dY t = (3 - 2uY t )dt + 2^Y t dW t and Y {k+1)St = f v^(l - wSt) + 



2St. 



The standard approach to analyze the time discretization error of the numerator and denominator of the 
left hand side of (|19p is then to introduce some functions v and w solutions to the partial differential 
equation: 



dtv = (3 - 2y)d y v + 2yd yy v - 6y\ (i, y) € R+ x (0, +oo) 



(20) 
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with initial conditions v(0,y) = y 2 and w(0,y) = 1. Now, we write (for the numerator, for example): 



E X|exp -9 / X*ds 



K 



x|.ex P -est^x 



k=l 




K-l 



fc-1 



J2 E («(T-fcJt,y fc «)-e- 



'u(T- exp -ft^Y" 



i-2 



fe=0 



J=0 



error bound of the form Cx5t can now be proved by some Taylor expansions as in [1, 16], provided the 
existence of a sufficiently smooth solution v to h20\) . We have not been able to prove existence of such a 
solution so far. 

1.2. Particle approximation 

We now introduce some notation to study the particle approximation. We recall that v denotes the 
number of large timesteps (the number of reconfiguration steps is v — 1), and At = nSt the time period 
between two reconfiguration steps. Let us suppose that we know the initial positions (X„ )i<j<jv of the 
./V walkers at time (n— l)At, for a time index n G {1, . . . , v}. The successive positions of the walkers over 
the time interval ((n — l)At, nAt) are then given by [X l n St , . . . , X l n KSt ), where (JQ t )o<t<At satisfies: 



X' 



x: 



, )0 + j\{X i n<s )ds+{w i 



t+(n-l)At 



w, 



(n-l)At 



(21) 



Here (W 1 , . . . , W N ) denotes a iV-dimensional Brownian motion independent from the initial positions of 
the walkers (X{ )i<i<N which are i.i.d. according to 2ipj(x)l{ x> o}dx. We recall that in our framework, it 
is possible to simulate exactly in law all these random variables (see Appendix) . We store the successive 
positions (X* st , , . . , X* nSt ) of the i-th walker over the time interval ((n — 1) At, nAt) in a so-called 
particle & G (K;) K (sec Figure [TJ: Vie {!,..., N},Vn G {!,..., v}, 



,x 



71. KSt 



(22) 



In the following, we will denote by £„ = . . .,£„ ) the configuration of the ensemble of particles at 
time index n. We have here described the mutation step. 




(n - l)At 



Figure 1. The i-th. particle £^ at time index n is composed of the successive positions 
(X l n st , . . . , X % n Kgt ) of the i-th walker on time interval ((n — l)At, nAt). 
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For a given configuration of the particles £ n at a time index n <G {1, . . . , v}, the selection step now 
consists in choosing the initial positions (X^ +1 )i<i<N of the N walkers at time nAt using one of the 
following resampling algorithm: 

(SI) The (X' l n+l )\<i<N arc conditionally independent w.r.t. and for 1 < i < N, X l n+l is dis- 
tributed according to the measure 

N 

^g(C n )S^ + (1 - e n g(C n )) d S d,^ (23) 
where g is defined by. for y = (y x , . . . , y K ) £ (K+) K , 

g(y) = cxp(-66tY / yi), (24) 



\ k=l / 

p> n denotes the weight of the j-th particle 

j 9 (Cn) fO^\ 

Ej=i .9(60 

and e n is a non negative function of such that e n < \j maxi<;<7v 9(£ n )- In particular the 
following choices are possible for e„: 

e n = 0, e„ = 1 or e„ = - — — . (26) 

maxi^^jv g{&) 

The so-called multinomial resampling method which corresponds to the choice e„ = gives 
rise to a maximum decorrclation with the former position of the particles, while with growing e n , 
more and more correlation is introduced. 

(52) The (Xn+i o)i<i<JV are such that 

\j g {i, . . . ,n}, v* g { (i + ECi 1 (Eii <4) } , 

X i = pi 

n+1,0 ^n 7 Ki (27) 

and the variables (X^ +1 o)i+yj™ x a < <i<N are conditionally independent w.r.t. £„, 
with o distributed according to Ejli {^Pn} K j [N - Eili a «) j 

where 

<= L^J. je{l,...,iV}. (28) 

Notice that the (X l n+l )i<i<N arc conditionally independent w.r.t. £„. This is the so-called 
residual resampling method. 

(53) The (^Q + i. )i<*<JV are such that, for 1 < % < N, 

N 

X (n+1),0 = P l n <{i-U* n )/N<Y. 3 l = 1 P l J^ n ' K ' ( 29 ^ 

i=i 

where (Ufyi<i<N are random variables i.i.d. according to the uniform law on [0, 1], independently 
of £„. Notice that the {X^ +1 o)i<i<JV are conditionally independent w.r.t. £„. This is the so-called 
stratified resampling method. 



For n € {1, . . . , u}, let us denote by 

1 N 



i=l 
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the particle approximation of the measure r\ n defined by: V/ : (R+) K — ► R bounded , 

E (/ (X(n-l)At+8t, ■ ■ ■ ) X(n-l)At+K8t) ex P (-^*Z)fel"l (XkSt) 4 

1»W = 7 7 in iw v\ — » (31) 



where the process (X t )o<t<T is defined by (|15p . 

For y = ( Vl , . . . , y K ) G (K;) K and / : (R* + ) K -> K, we set 

P/(y)=E(/(^,...,^ t )) (32) 

where for x G R?j_, 

X? = x+ [ b(X*)ds + W t (33) 



o 



denotes the solution of the stochastic differential equation (fl"5)) starting from x. By the Markov property, 
the measures (?7n)i<n<y satisfy the inductive relations, for any function / : (R^) fi — > R bounded, Vn € 



Vn+l(f) 



E exp (-^Efc=i(^«) 4 ) E / (X„ A t+«, ■ ■ -,W. 



>c<5i ) 



P^<5t) 



(34) 



J7 n ( ff )E (exp (-^Efe="i 1)K (^* 
1 E (gP/ (X (n _ 1)At+5t , . . .,X {n _ 1)At+KSt ) exp (-^EiTi 1)K (^t) 4 )) T) n (gPf) 



^9) E(exp(-^Ei="i 1)K (^t) 4 )) ^ 

(35) 

where g is defined by (|2~4")) . Moreover, we can express E^ MC (T) defined by (fTB]) as: 

J# M c(T) = §« + *^^. (36) 
2 IMS') 

Therefore the particle approximation of Eumc {T) is given by 

eZc K (T) = ^ + 9 V 1^. (37) 

This approximation of -Bdmc (T) corresponds to the expression (JSj) given in the introduction. We will also 
prove in Corollary[13]below the convergence of the approximation which corresponds to the expression (|7|) 
given in the introduction (see Equation (|4T))) below). 

The convergence of the approximation E^-^£(T) is ensured by our main result : 
Theorem 4. 

E 



E DMC (T) - EZ'ciT) < — + (38) 



C_ + Cu 

UK y/N' 

where the constant C only depends on T and the constant C v on T and v. 

Remark 5. The number of selection steps is v — 1. For instance, when v = 1, there is no selection 
involved in the expression of (T) and the particles remain independent. In this case, the first term 

in the right hand side of i38\) corresponds to the time discretization error proved in Proposition QJ while 
the second term is the classical error estimate related to the law of large numbers. For a fixed number 
of selection steps, the theorem ensures the convergence of the particle approximation E^^£(T) as the 
time-step St = T / (vk) used for the discretization of the stochastic differential equation (|15l) tends to 
while the number N of particles tends to +00. But this result does not specify the dependence of C v on v 
and gives no hint on the optimal choice of the number of selection steps in terms of error minimization. 
We are going to deal with this important issue in the numerical study (see Section^). 
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According to the above expressions ([36]) and fl3Z]) of E$ MC (T) and E^'^(T), this theorem is easily 
proved by combining Proposition [1] and the following result : 



Proposition 6. 



E 



Proof of Proposition [6] : One has 



€(9) vAg) 



< ■% (39) 



E 



' < n I! ?(gyt)-vM)\ 



2 



fiAg) 



?(gyt)V\ 1/2 ( E (v?(g)-v»(g)r 



According to Proposition [7] and Lemma [T2l below, the first term of the right-hand-side and the quotient 
in the second term are smaller than C„/y/N. Since by Jensen's inequality, ( j l '^ l N 9 ^ S j < ~^fsF'' ^ e 
boundedncss of E ( j^jw^^- ) follows from Lemma [5] below. I 



Proposition 7. For any bounded function f : (R!j_) K — > BL 

Vn G {1, . . . , i/}, E(fa*(/) - Vn (f)) 2 ) < (40) 
where the constant C n does not depend on k. 

For any function f : (R^) K — > R swc/i i/ia£ for some p > 2, ||/|| K . P = sup ^^j, is finite, 

Vn G {1, . . . , u}, E\ v Z(f) r ln (f)\ < %11/lkp, («) 



where the constant C n does not depend on k. 

For / bounded, the first estimate ([40|) is proved in [9]. In order to prove Proposition O we need to 
apply Proposition [7] with f(y) = g(y) and f(y) = g(y)y^, which are bounded functions with L°° norm 
respectively equal to 1 and ^ where C is a constant not depending on 5t. But we want to obtain 
the convergence when St tends to 0. This is why we need the second estimate (|4Tj) . that we use with 
f(y) = g(y)yf. for which ||/|| K ,p is bounded and does not depend on St. 



Notice that for / bounded, Corollary 2.20 in [9] states the convergence in law of s/N(r]^ (/) — Vn(f)) to a 
centered Gaussian variable and gives an expression of the variance of this limit variable. Because of the 
complexity of this expression, using this result with f(y) = g(y)y^ did not really help us to understand 
the dependence of C v on v (sec Remark [5] above) . 

Proof : For / bounded, the first estimate (|40|) is proved by induction on n in [9] (see Proposition 2.9). 
Since we follow the same inductive reasoning to deal with / such that ||/|| K . P < +00, we give at the same 
time the proof for / bounded. 

Since the initial positions (£i)i<i<jv are independent and identically distributed with £\ K distributed 
according to 2i(jj(x)l{ x>0 }dx, the statement holds for n = 1. 

To deduce the statement at rank n + 1 from the statement at rank n, wc remark that according to (|35[) . 

«^+i(/) - Vn + i(f) = T n+1 + -L^ ((v%(gPf) - Vn(gPf)) + ^y/W g) - € (g))) (42) 
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where we recall that P is defined by (|32|) . and 



T n +i = Vn+i(f) - 



To deal with this term T„+i, one remarks that for the first type of selection step (SI), all the possible 
choices of e n given in (p?!)]) are <7(£ n )-measurable. As a consequence, for i £ {1, . . . , N}, 

N 

E(/(£+i)|£„) = en9(C)Pf(&) + (1 - e n g(£))'£fiPf(£), 

3=1 

where is defined by (|25p . Multiplying this equality by and summing over i, one deduces 
Now, for the stochastic remainder resampling algorithm (S2), by (|27| . E(r?^j_ 1 (/)|^ n ) is equal to 



(43) 



1 W 

3=1 



^(a 



p/(&) + E 



1 < ^ ^ 



and (|43|) still holds. Finally, for the stratified resampling method (S3), by (|29|) . we have (using the 
footnote 2 ) 



AT JV 



t=l 3=1 

1 w / w 

3 = 1 \i=l 



l <(i-^A)/Jv<EL 1 pi,} 



1 N 
-y , 



3 = 1 



;=i 



3-1 

^E^ + c/ « 

i=i 



Cn) Pf{&), 

tn) p/(e), 

6, ) P/(£), 



= I>„p/(&), 

3=1 

which yields again (|4"3"| . Since for all three possible selection steps, the variables (£^ +1 )i<i<iv are inde- 
pendent conditionally on £„, one deduces that 

1 N 1 

E((T Il+1 ) 2 ie„) = t^E e - E(/(c\ + i)iCn)) 2 leJ < ^e . 



i = l 



Therefore 



E((T„ +1 ) 2 )<-E(^ +1 (/ 2 )). 



(44) 



When / is bounded, v^+iiP) < 11/11 



'),, (gP/) 



< ||P/|U and ||P/|U < H/lloo. Hence by 62]), 



E((^ + i(/)-^+i(/)r)<3(^f i + 



/HSo , E((^(.gP/)- ? y„( g P/)) 2 ) + ||/|| 2 E((^(.g)-, 7 „( 5 )) 



(^(ff)) 2 



with the second term of the right-hand-side smaller than CWfW^/N by the induction hypothesis and 
Lemma [12] below. 
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Kl n m m , ^(« +1 (/ 2 ))) 1/2 , El^^P/)-^^/)! 

E - *7n+l(/)| <" 



3 



n»(j>Pf) \ 2 \ 1/2 W(g)-^n(g)) 2 ) 1/2 



Since ||/ 2 || fe , 2p < 2||/|| 2 p (by using the inequality f 2 {y) < 2\\f\\ 2 K _ p (l + the first term of the right- 

hand-side is smaller than C n \\f\\ KlP /y/N by Lemma [9] below. Since, according to Lemma [TOl below. 
II-P/II k,p < e c ' ! ' At ||/|| Ki p, the second term is smaller than Cnll/IUjp/V'N by the induction hypothesis and 
Lemma 1121 Last, by using successively Cauchy Schwartz inequalities, (|43[) for f 2 and Lemma [5J one 



obtains that E (^gf^) 2 < E ( ^fgfi ) < E (^^P) = «+i(/ 2 )) < ^ 
Proposi 

c n /V]v. 



from the Proposition statement for / bounded and Lemma [12] that ( ^L_?^£^ ) — j s smaller than 



Remark 8. Proposition^ (and therefore Theorem^ also hold for the stratified remainder resampling 
algorithm, which consists in combining the stochastic remainder resampling and the stratified resampling. 
More precisely, it consists in replicating \_Np l n \ times the i-th particle, and then completing the set of 
particles by using stratified resampling to draw the N R = N — Ei=iL-^/°nJ remaining particles, the i-th 
particle being assigned the weight p R,z = {Np' l n }/N R . 



Lemma 9. Let h : (R* + ) K — ► R + be such that for some p > 2, \\h\\ KtP < +00. Then, 



<^" At ii^iu,p(i+E(x n, 



where Xq is distributed according to the measure 2ip 2 (x)l{ x> o}dx (see 115\) ). 



Proof : As the variables K , 1 < i < N are distributed according to the invariant measure 2ip](x)l{ x>0 }dx, 
one has E(?yf(/i)) < ||ft|U, p (l +E(X ) P ). In addition for n > 1, according to ((43]), E(r)% +1 (h)) = 
E ( j'^Nfefi j where ||Pft.|| KiP < e c ' J,At |j/i||fc j p by Lemma [TUl Therefore it is enough to check the bound for 

E 



(gh) \ 



(g) J' 
For n > 0, one has 



[^^-j<\\h\\ K<p \l + E\ — V ■ _ . ) _ I 



Ef =1 exp(-^ELi(^ +1 , fc ) 4 
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Let us denote in this proof = ^Q+i oj where < n < v — 1 and 1 < i < N. Let us set 
F = fc> 1 < « < < k < k — 1). By Lemma E] below, 



El ~ 



^ < 



V Ef=iCx P (-^ELi(^ +1 , fc ) 4 ) / " Ef = iCxp(-^Efc=!(d, fc ) 



E*=i exp ( -^E^iCC+i,*) 4 ) E((Xf t )f)| x= ^ i 



Ef=iexp(-^E^(e +1 



)4 



r ^Etiexpf-^E^^+i,,.) 4 ) (C+i,k-i) p c 



where we have used the definition of the mutation step (see ([21])) and the Markov property for the stochas- 
tic differential equation ([55)) to obtain the equality, and then Lemma[TU]for the last inequality. Notice that 

this estimate also holds for n = 1, in which case the right hand side reduces to (Cn+i o) P + e Cpdt — 1. 
Taking expectations and iterating the reasoning, one deduces that 

/E 4 =iCxp(-^ELi(a +1 , fc ) 4 )(C; + i. K ) p \ e c P A t J^ S^r™ 
V Ef=iexp(-^ELi(e +M ) 4 ) / * ti to 



Inserting this bound in (|4"5)) , one concludes that 



< e u ^\\h\\ K , p 1 + E 



jvE^+i.o) P jl • 



For n = 0, one deduces that E ( ~tt^ ) _■ e CpAt ||/i|| K . p (l +E(Xq)), where Xo is distributed according to 
the measure 2?/>f (a^l^x^dcc. 

»7n {g{y)yl) 



( l N \ (n N i 

For n > 1, since by a reasoning similar to the one made to obtain (|4"5|) . E I — / J(£n+i o) P I = ^ ( — ^ 

V i=i / ^ ?? « 



one also deduces that 



E (t7^) < e u "^\\h\\ K , p E 
The proof is completed by an obvious inductive reasoning. 



Lemma 10. For any p>2, there is a constant C p such that 

Vx g Vt > 0, E((Xf ) p ) < (1 + x p )e Cpt - 1, 

where Xf is defined by H33\) . Therefore, if h : (M!j_) K — > K is swc/i i/iai ||/i||«, p < +oo then \\Ph\\ KtP < 
e C,,At||^|| where the operator P is defined by 



Proof : By Ito's formula, d{Xf ) p = (*£^±_) (Xf ) p ~ 2 - wp(X?)A dt + p(XfY~ 1 dW t . Hence 

(x?T < * + jf + p(p + 1 2 ' 2a;) (^T) ds +P j\x°)^dW s . 
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Formally, taking expectations in this inequality, one obtains 

K(x ( ?) < ^ + 1 ^ + p{p+ \- 2uj) n(x*T)d S , 

and check by Gronwall's lemma that the conclusion holds with C p = p £ • This formal argument can 
be made rigorous by a standard localization procedure. 
For h : R* — > R such that < +00 one deduces that 

Vy 6 R+, \Ph(y)\ < E\h(X y s ?,. . . , X^)\ < C\\h\\ K , p (l + E((J^)) < e c > M \\h\\ K , p (l + £). 



Lemma 11. 

V(zi,... ) * w ),(oi,...,Oi V )GR^ rift ^a,>0, Vp>0, Vc > 0, ^ 4 ' < 



Proof : Let us set /(c) = i=1 ' ' _ — j-. By Holder's inequality, the derivative 



/'(c) = 



i=i Oi«i e cz * Ei=i <»i«i e cz ' \ Ei=i Oi< e 



EAT _„4 v-^JV _„4 I v — i/V 



is non positive. Hence for any c > 0, /(c) < /(0) = ~i.lv °' Zi . 

Lemma 12. 77ie sequence (j]n{g))i<n<u is bounded from below by a positive constant non depending 



Proof : Since 

Cg) = E(exp(-^E^ 1 ^ t )) 

E^xp^Eir! 1 ^^" 

the sequence (??n(<7))i<n<i/ is bounded from below by 



l[r 1n (g)=Eiexpl-65tJ2Xk S t 

n=l \ \ k=l , 



According to Lemma [2J this expectation converges to E ^cxp \ —9 j Q Xjdsjj > when k tends to +00, 
which concludes the proof. 



We can now prove, as a corollary of Theorem [4] the convergence of the approximation E D -^£(T) of 
Ebmc{T), defined by: 

3 



N 



dmc 1 



-OH 

2 N 



(46) 



i=i 



Corollary 13. 



E 



Eumc{T) - £dmc( t ) 



< 



C C v 



where the constant C only depends on T and the constant C v on T and v. 
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Proof : By using the result of Theorem |4] and Cauchy Schwartz inequality, it is sufficient to prove the 

\ 2 , c v 



estimate E ( E^-^£(T) — E£^q(T) ) < Let us denote in this proof = X % u+l for 1 < i < N. 
We have: 



N 



N 



DMC K 1 ) DMC V J I 



/ 1 N \ 1 N 

\ i=l J i=l 



by using the fact that, for any function / 



/ 1 N 



€{9(v)f(v,)) 



(47) 



which is obtained by a reasoning similar to the one made to prove (|43[) . Now, using the same method as 
to obtain (|4"4")) . one easily gets the estimate: 



E (*£&"(T) - E^{T)) < °-E l± £(£ +1 , ) 8 J = 



by using again (|4T|) . Lemma [5] completes the proof. 



We end this Section by proving that Proposition [6] also holds for the numerical scheme (fl7|) . 

Proposition 14. Lei us consider the Markov chain {Xjgt)o<j<K generated by the explicit scheme (|17j) 
and denote by Q its transition kernel. We now define the measure rj n by replacing (Xjst)o<j<K with 
(Xjst)o<j<K in 131]) . and we define accordingly the evolution of the particle system: conditionally on £ n , 
the vectors (X^ +10 , X^ +1 st , . . . , X^ +1 K st)i<i<N are independent, with K+i o)kkjv distributed accord- 
ing to the selection algorithm (SI) (see ((23)) ). (S2) (see (|27|) ) or (S3) (see (|29|) ). and (X^ +1 jg t )o<j<K a 
Markov chain with transition kernel Q. Then, we have: 



3 



vi!{gyi) vAgyi) 



< 



c v 



N 



Proof : Looking carefully at the proof of Proposition [5] above, one remarks that (f3T)|) holds in this 
framework as soon as Lemma [T"2l holds, and the following property, which replaces Lemma [TOl is satisfied: 



3C > 0, Mx e R+, Qf(x) < e cst (l + f(x)) - 1 for f(x) = a; 4 and f(x) = x* 



(48) 



2 1 /2 

Let us first prove (gSJ). We have: Q/(s) = E (/ (Xf t )) where Xf t = ((1 - Lu5t) 2 x 2 + 2xW St + (1 _3 t)2 + 2<ft) . 
Now, for q G N* , 



(X£) 2 « = ^ 



Ji!j2!j3! 



(l-£j5t) 2jl 2 J2 x 2j ' 1+J2 W\ t 



2St 
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where the indices (ji,j2, J3) are non negative integers. Remarking that the expectation of the terms with 
32 odd vanishes and then using Young's inequality, one deduces that for 8t < w-, 



E((X^) < (1-c^)V+e(^- 



2St 



- UjSt) 2 

< x 2 " + C q 6t + C q (x^'St + St^i 1 - 

31+32+33=9 

ji<q,j2 even t js<q 

< (1 + C q 8t)x 2q + C q St < e c « St (l + x 2q ) - 1. 



C q J2 * 2(9 -™ } <5 



31+32+33 = 9 

ji<1,j2 even ,j 3 <q 



32 + 233 , 



(49) 



Let us now prove LemmalT^lfor the scheme p7|) . As noticed in the proof of LemmafTZlabove. it is sufficient 
to bound from below E (exp \—65t Y^k~i ^tst)) ■ By Jensen inequality, we have E (exp {—68t Y^k=i ^kSt)) — 
exp (-0^ 22,1=1 E ( X kSt))- B y usin S ©i {t is eas y to P rovc b y induction that E (X% st ) < e C2kSt {l + 
E (Xq)) — 1 and this concludes the proof of Lemma [T2l in this framework. 



In order to obtain a complete convergence result of the form (f38|) for the scheme p7|) . it remains to prove 
the complementary bound (|19p . that we have not obtained so far. However, we will check by numerical 
simulations that Q38|) still holds. 



2. Numerical results 



2.1. Computation of a reference solution by a spectral method 

In this section, we would like to explain how we can obtain a very precise reference solution by using a 
partial differential equation approach to compute -Edmc(^) (see [3]). 



2.1.1. A partial differential equation approach to compute Edmc(T) 

Let us introduce the solution to the following partial differential equation for : 

I ^ = -H<t>, (t,x)eR+xR (50) 
\ </>(0, a;) = ipi(x), i6l 

where H (resp. ijjj) is defined by (fTD]) (resp. (HH)). Since ipj £ H, it is a standard result that this problem 
admits a unique solution G C°{R + ,H) n C°(R* + , D H (H)) n C 1 (K+,H)- The function is regular and 
odd, and therefore is such that <fi(t,0) = for all t > 0. Therefore the function is also solution to the 
following partial differential equation: 

^ = -H6, (t,x) £l + xl 
(j)(t,0) = 0, t > 
0(0, ar) = ^/(x), i£l 

In [3], we have shown that since satisfies (f5"Tj) . we can express -Edmc(0 (defined by using the 
function (see Proposition 11 in [3]): 

_ (ff^,0(f)) 

<W,0(*)) 



Our reference solution Bdmc(^) w iU rely on formula (|52|) after discretization of (|50"|) by a spectral method. 
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2.1.2. Computation of the wave function <f> 

We will briefly present the spectral method developed to compute an approximation of 4>. We recall that 
the Hcrmitc polynomials are defined by : 

VneN, h n (x) = {-l) n e x2 ^{e-* 2 ). 

We introduce the eigenfunctions of the operator Hq, normalized for the L 2 (M.) norm associated with the 
eigenvalues E n = ui(n + 1/2) for n > 0, 



(p n (x) = h n (*Suix) exp(-iwx 2 ) y 



(^A) 1/4 



It is well known that the vector space spanned by the set of functions {f2k+i}k>o is dense in Vo = {<p G 
H 1 (M.) n H | xcp G L 2 }, which is the domain of the quadratic form associated with Hq. 

Let us now introduce the functional space V = {p E H 1 (R) n Ti \ x 2 Lp £ L 2 }, which is the domain of the 
quadratic form associated with H. The set of functions {tp2k+i}k>o is also a basis of V. 

Let V n — Span((pi,ip3, . . . ,ip2n—i)- We use this approximation space to build the following Galcrkin 
scheme for (|5U]t : find <f> n G C°(E+, V„) such tha10 </>„(0, x) =ipi, and VV G V„ 

^P-,<p\=-(Hct> n (x,t),<p). (53) 

We diagonalize the operator iJ restricted to V n . We denote (yoiVai • ■ ■ i Vn-i) the eigenfunctions and 
£^0; £"2 j ■ ■ ■ i ^n-i ^ ne associated eigenvalues. Because of the symmetry of H, it is easy to check that V n 
can also be spanned by ((^q j V2j ■ ■ ■ > <Pn-i) : 

V n = Spanitf, ft^). (54) 

Since for t > 0, n (t, .) G V ra , there exists Mfc(t), k = 0, . . . , n — 1, such that 



$>*(*K" (55) 



fc=0 



In view of (|54|) and (|55|) . (|53[) is equivalent to the equations: Vi = 0, . . . , n — 1, 

fc=0 \ fc=0 / 

n-1 



fe=0 



Wc deduce that V7c = 0, . . . ,n — 1, 



so that 



where u fc (0) = (-0/, 



b n (t,x) = Y / M0)exp(-E^l(x), (56) 



fc=0 



3 Notice that ifij = <p\ e V„ 
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Remark 15. The eigenf unctions of H are obtained by diagonalization of the matrix A = (aij )i,j=o,...,ra— 1 
with Vi, j = 0, . . . , n — 1 : 

= (H(f2i+i, <P2j+i) , 

= (H Q Lp 2l+ l,ip2j+l) + d (x i ip 2i+ l,ip2j+l) , 

= 5ij w (2i + §) + 8 (x 4 (p 2i+ i,(p2j+i) ■ 

We can use the n-point Gauss-Hermite formula to deal with the integration of the second term on the 
right. We recall that this method provides an exact result for f_°° p(x) exp(— x 2 )dx as long as p is a 
polynomial of degree 2n — 1 or less. 

2.1.3. Approximation of Eumc(T) 

We now use formula (|52[) to approximate Edmc(T). By an elementary calculation, we obtain the following 
approximation: 

Edmc (T) * . ( 5 7) 

In our test cases, we have observed that n — 40 is enough to reach convergence. 

Notice that for a given n, the convergence in time to the lowest eigenvalue Eq is exponentially fast, with 
an exponent equal to the spectral gap — Eg. 



2.2. Numerical results of Monte Carlo simulations 



In this section, we perform various numerical experiments to validate our theoretical results, and to 
explore some features of DMC computation. In particular, we propose in Section 12.2.21 an empirical 
method to determine the optimal number of reconfigurations. In all the computations, the final time is 
T = 5, which appears to be sufficiently large for the convergence t — > oo to be achieved with enough 
accuracy. 



2.2.1. Error and variance as a function of the numerical parameters 

We represent on Figured the expectation e and the variance v of the error : -E^dmc * CO — ^dmc (T) 
as a function of the number of walkers N, the time step St and the number of reconfigurations v — 1 
where Edmc(T) is approximated using (j57]l and E^^ /{vSt \T) is defined by fl£fl. The multinomial 
resampling method (which is (SI) with e„ = 0) was used. 

The top figures represent the expectation of the error and its variance according to the number of 
walkers. To compute these quantities, we perform 2000 independent realizations, with the number of 
reconfigurations v — 1 = 50, a small time step St = 5.10 -3 and 9 = 0.5. The simulations confirm the 
theoretical result : the error decreases as C/\N. 

The effect of the time step is shown on the two figures in the center. The numerical parameters are: a 
large number of particles ./V = 5000, number of configurations v — 1 = 30, 6 = 2 and 300 independent 
realizations. We can see on the figure on the left that the error decreases linearly as the time step 
decreases. We also remark that the error is smaller with the approximate scheme f| 1 T[) than when using 
the exact simulation of the SDE (|15p proposed in the Appendix. This rather amazing result can be 
interpreted as follows. When using the exact simulation of the SDE, there is only one source of error 
related to the time discretization, namely the approximation of the integral in the exponential factor 
in P]). When using the scheme (|17p . we add a weak error term which seems to partly compensate the 
previous one. 




Figure 2. Expectation and variance of the error when (fT5|> is discretized according to 
the method described in Appendix (dotted curve) and according to the scheme (fT7|) (solid 
curve) . 



The last figures represent the effect of the number of reconfiguration steps. The numerical parameters 
are: time step St = 5.10 -3 , number of particles N = 5000, (9 = 2 and 300 independent realizations. The 
curve representing the variation of the error according to the number of reconfigurations has the shape of 
a basin. We deduce that on the one hand a small number of reconfigurations has the disadvantage that 
walkers with increasingly differing weights are kept. On the other hand a large number of reconfigurations 
introduces much noise. An optimal number of reconfiguration seems to lie between 20 and 50. 

2.2.2. Optimal number of reconfigurations 

On Figure [3J we check that the optimal number of reconfigurations in terms of the variance v of 
S DMC /( " ft) ( T ) ( and not of the error as in Section 12.2. lj) is also obtained for a number of reconfigu- 
ration which seems to lie between 20 and 50 (using again the multinomial resampling method). The 
numerical parameters are those considered for the figures below in Figure [2j time step St = 5.10~ 3 , 
number of particles N = 5000, = 2 and 300 independent realizations. We have not studied how the 
optimal number of reconfigurations varies according to the other numerical parameters. 

We have investigated a practical method to estimate numerically the optimal number of reconfigurations. 
On Figure|3]we represent the variance of Edmc^'W according to time t, without any reconfiguration step 
(which corresponds to v = 1). The other numerical parameters are again those considered for the figures 
below in Figure [2j We observe that the variance is minimal at t* ~ 0.25. We remark that v = T/t* = 20 



TITLE WILL BE SET BY THE PUBLISHER 



21 




Figure 3. Variance of E D ^ C (T) in function of the number of reconfigurations 
when (|15p is discretized according to the method described in Appendix (solid curve) 
and according to the scheme (fTT)) (dashed curve). 



is close to the optimal number of reconfigurations obtained on the previous figures. We have checked this 
empirical result for various sets of the parameters. It seems that the optimal number of reconfigurations 
is related to T/t* where t* minimizes the variance of E^^/ St (t). Since v = 1, no selection step occurs 
and the particles are thus independent. According to the multidimensional central limit theorem, the 
variance of E^^'^ St (t) can be approximated by 

1 / Var(F t ) Covar(y t ,Z t ) a Vax(Z t ) 

N {WW " { t) WZtff + { { t)) (W 



where 



and 



t/st 

Y t = E L (X t )exp I -5tY^E L {X kS t) 



fc=l 



t/St 

Z t =cxp I -5tJ2E L (X kSt ) 



fe=i 



Therefore, the optimal number of reconfiguration steps could be estimated by this method, through a 
precomputation over a few independent trajectories. 



2.2.3. Comparison of the resampling algorithms 

We finally compare various resampling algorithms on Figure El where the variance of E^-^'^^ St \t) as a 
function of time is represented. The numerical parameters are: N = 1000, St = 5.10~ 3 , v — 1 = 20, 9 = 2 
and 200 independent realizations. 

We first observe on the figure on the left that without any resampling, the variance of the results explodes 
with increasing time. This shows the necessity to use resampling algorithms. We compare the follow- 
ing resampling algorithms: multinomial resampling (which is (SI) with e„ = 0), correlated multinomial 
resampling (which is (SI) with e„ = 1/ maxi<.;<jv <?(£«)), residual resampling (which is (S2)), stratified 
resampling (which is (S3)), stratified remainder resampling (which combines residual and stratified re- 
sampling, see Remark |8]) and systematic resampling (which corresponds to stratified resampling with 
U} % = ... = = U n , see the Introduction). We observe that, as expected, when more correlation 
is introduced, the variance due to the resampling is reduced. The multinomial resampling method is 
generally the worse, while the best resampling methods seem to be systematic resampling or stratified 
remainder resampling. 
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Figure 5. Variance of £^mc (*) 

as a function of time t, for various resampling 
algorithms: Without = without resampling, Mult = multinomial resampling, CMult = 
correlated multinomial resampling, Res = residual resampling, Strat = stratified resam- 
pling, StratRcm = stratified remainder resampling, Syst = systematic resampling. 

Conclusion 

In this paper, we have proved on a simple example convergence of numerical implementations of the 
DMC method with a fixed number of walkers. The theoretical rates of convergence are confirmed by 
numerical experiments and are likely to hold in more general situations. We have also checked numerically 
the existence of an optimal number of reconfiguration steps. Various resampling algorithms have been 
considered, both theoretically and numerically. For future work, we plan to investigate criteria devoted 
to the choice of the number of reconfiguration steps. One interesting direction is the use of automatic 
criteria based on a measure of the discrepancy between the weights carried by the walkers to decide when 
to perform a reconfiguration step. 



Appendix : Simulation of the stochastic differential equation (fl5j) 

In this appendix, we show that it is possible to simulate exactly in law the (if+l)-plet (^0; X$t, ■ ■ . , Xxst), 
where X t is defined by ()15p . Let (G, U) denote a couple of independent random variables with G normal 
and U uniformly distributed on the interval [0, 1]. 
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Simulation of the increment X t — X s , for t > s. 

The square Rt of the norm of a 3-dimensional Brownian motion Wt = (Wj, W£,Wf) solves dR t 

/•* w s ■ dW s R t 
3dt + 2\/R~ t dB t where B t = I — | s - is a one-dimensional Brownian motion. Hence p t 



,\W S \\ ' l + 2ut 

solves 

dp t = {3-^)^+2^-^. (58) 

It is easy to check that ^ J 2 " ^ ^ -^=t=^ is a Brownian motion. Hence, performing a time-change 
in (|58|) . one obtains that pj_( e 2ut_i) = e~ 2u ' t Rj_( e 2u,t_ 1 } is a weak solution of the equation dYt — 
(3 - 2uY t )dt + 2^JY t dW t satisfied by Y t = Xf. Therefore e~ ut . /r~j~Z^~^ is a weak solution of JTSJ). 



j(^ 2 

For u > u, R v has the same distribution as (y/EZ + W\ - W^f + {W 2 V - W 2 ) 2 + (W 3 V - Wf t ) 2 , and 
therefore as (\/R^+ G\/v — u) 2 — 2(v — u) \og(U) with (G, U) independent from R u . Hence for t > s, X t 
has the same distribution as 



2 \ \ 1 / 2 

e us I s + -!^=(e 2wt - e 2 ^ 8 ) 1 ' 2 ] - 2 — (e 2uJt - e 2uJS ) \og(U) ) ) 
\/2uj J 2lo 



where the couple (G,U) is independent from JT S . 



1/2 



Simulation of Xq with distribution 2tpj(x)lr x>0 xdx. 

1/2 

The random variable -j== {G 2 — 2 log(?7)) is distributed according to the invariant measure 2ipj(x)l{ x>0 }dx, 
as suggested by letting the time increment t — s tend to +oo in the previous simulation. Indeed, 
G 2 - 21og(E7) is a Gamma random variable with density 2 3 / 2 r(3/2) ^{z>o}\^ e ■ And one deduces 

1 /2 

the density of -k= (G 2 — 21og(t/)) by an easy change of variables. 
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